/*----defines the maximum size of K*K matrix.
	This parameter and these objects involve in "MultiplySquareMatrices" and "RWishart"-----------------------------------------*/

#define	Klimit	50
#define Pi 3.14159265358979323846
double	Intermediate[Klimit*Klimit];
double	Lowertri[Klimit*Klimit], Vi[Klimit*Klimit], sqrtVi[Klimit*Klimit], Nij[Klimit*Klimit];


/*---- define structures --------------------------------------------------------------------------------------------------------*/
/* Structure of Y */
typedef struct
{
	double	*observations;			/* observations */
	double	*errors;				/* residual errors */
	double	*stobs;					/* standardized observations */
	double	*meanerrors;			/* mean residual erros (for DIC calculation) */
	int		nmiss;					/* number of missing records*/
	int		*whichrec;				/* which records are missed*/

}	Ystruct;

/* Structure of location parameter */
typedef struct 
{
	double	*covariates;			/* values of covariates */
	double	x2;						/* sum of squared covariates */
	double	currenteffect;			/* current effect */
	double	*sampledeffect;			/* sampled effect */
	int		nrec;					/* number of non-missing records that have the effect */
	int		*whichrec;				/* indicate which record has the effect */

}	Lstruct;

/* Structure of Lambda */
typedef struct  
{
	int		update;					/* model for Lambda */
	int		down;					/* downstream (affected) trait */
	int		up;						/* upstream (affecting) trait */
	double	*lambda;				/* lambda value for each record */ 
	double	*sampledlambda;			/* sampled lambda */
	double	meanlambda;				/* mean of the prior distributio of lambda */
	double	lambdasigma;			/* variance of the prior distribution of lambda */
	double	*lambdahat;				/* mean of the posterior normal distribution */
	double	*vlm;					/* variance of the posterior normal distribution */
	double	*y2;					/* y^2. When model 1, sum of all y^2, when model 2, y^2 of each record*/
	int		*lambdaclass;			/* classes of records, only when 2*/	
	int		*nic;					/* number of records included in each class, only when model 2 */
	int		nc;						/* number of classes, only when model 2*/
	double	*lambdalist;			/* lambda value unique to each class, only when model 2*/
	double	alpha;					/* concentration parameter, only when model 2 */
	int		*sampledclass;			/* sampled classes, only when model 2 */
	int		Nsp;					/* Number of splines. Only when model 3*/
	int		Ord;					/* Order of the base function. Only when model 3 */
	int		*HaveithSP;				/* Records with ith spline > 0 */
	int		*NhaveithSP;			/* Number of records with ith spline > 0 */
	int		Naccepted;				/* Number of acception of MH update. Only when model 3 */
	double	*B;						/* Base function. only when model 3 */
	double	*P;						/* Spline coefficients. only when model 3*/
	double	*sampledPi;				/* Sampled coefficients. only when model 3*/
	double	*knot;					/* Positions of the knots */
	double	a;						/* Interval of the knots. only when model 3 */
	double	sigmaPi;				/* Prior variance of the spline coefficients */
	double	Nupi;					/* A hyperparameter of the prior inverse chi square distribution of sigmaPi */
	double	S2pi;					/* A hyperparameter of the prior inverse chi square distribution of sigmaPi */

}	Lambdastruct;


/*----For random number generator------------------------------------------------------------------------------------------------*/
int ix, iy, iz;

/*----initialize the random seed (by STRUCTURE)----------------------------------------------------------------------------------*/
int Randomize ( int i )
{
	if ( i == -1 )
		i = (int) time ( NULL );
	srand ( i );
	return ( i );
}

/*----generate a random real number between 0 and 1 (by Ayres1998)---------------------------------------------------------------*/
double WichHill()
{
	ix = (171*ix)%30269;
	iy = (172*iy)%30307;
	iz = (170*iz)%30323;
	return ( fmod ( ( ix/30269.0 ) + ( iy/30307.0) + ( iz/30323.0), 1.0 ) );
}

/*----generate a random real number ---------------------------------------------------------------------------------------------*/
double  RandomReal ( double low, double high)
{
	return ( low + WichHill() * (high - low) );
}

/*----generate a real number from Uniform (0,1) (by STRUCTURE with some modifications)-------------------------------------------*/
double rnd()
{
  double value;

  do
    value = WichHill();
  while ((value==0.0)||(value==1.0));

  return value;
}

/*----generate a real number from Gamma (n, lambda) (by STRUCTURE), E[x] = n/lambda, V[x] = n/lambda^2 --------------------------*/
double RGamma(double n,double lambda)
{
  double ar;
  double w;

	double x = 0.0;
	if ( n < 1 )
	{
		const double E = 2.71828182;
		const double b = ( n + E ) / E;
		double p = 0.0;
		one: 
		p = b * rnd ();
		if ( p > 1 ) goto two;
		x = exp ( log ( p ) / n );
		if ( x > -log ( rnd() ) ) goto one;
		goto three;
		two: 
		x = -log ( ( b - p ) / n );
		if ( ( ( n - 1 ) * log ( x ) ) < log ( rnd() ) ) goto one;
		three:;	
	}
	else 
	{
	    if ( n == 1.0 )
		{                         /* exponential random variable, from Ripley, 1987, P230  */
		    double a = 0.0;
		    double u, u0, ustar;
            ten:
		    u = rnd();
	     	u0 = u;
         	twenty:
		    ustar = rnd();
		    if ( u < ustar) goto thirty;
	    	u = rnd();
		    if ( u < ustar) goto twenty;
	    	a++;
		    goto ten;
        	thirty:
		    return ( a + u0 ) / lambda;
		}
	    else
		{
	    	double static nprev = 0.0;
	     	double static c1 = 0.0;
	    	double static c2 = 0.0;
	    	double static c3 = 0.0;
	    	double static c4 = 0.0;
	    	double static c5 = 0.0;
	    	double u1;
	    	double u2;
	    	if(n != nprev)
			{
		    	c1 = n - 1.0;
		    	ar = 1.0 / c1;
		    	c2 = ar * ( n - 1 / ( 6 * n ) );
		    	c3 = 2 * ar;
		    	c4 = c3 + 2;
		    	if( n > 2.5 ) c5 = 1 / sqrt(n);
			}
		    four:
	       	u1 = rnd();
	       	u2 = rnd();
	       	if ( n <= 2.5 ) goto five;
	       	u1 = u2 + c5 * ( 1 - 1.86 * u1);
	       	if ( ( u1 <= 0 ) || ( u1 >= 1 ) ) goto four;
	       	five:
	       	w = c2 * u2 / u1;
	        if( c3 * u1 + w + 1.0 /w < c4) goto six;
		    if( c3 * log ( u1 ) - log ( w ) + w >= 1) goto four;
	       	six:
	       	x = c1 * w;		
	       	nprev = n;
		}
	}

	return x/lambda;
}

/*----generate a real number from chi square distribution whose degree of freedom = def -----------------------------------------*/
double Rchi( double def )
{
	return ( 2.0 * RGamma ( def / 2.0, 1.0 ) );
}

/*----generate real numbers from N (0,1) ( by STRUCTURE )------------------------------------------------------------------------*/
double snorm()    /*was snorm(void) -- JKP*/
{
static double a[32] = {
    0.0,3.917609E-2,7.841241E-2,0.11777,0.1573107,0.1970991,0.2372021,0.2776904,
    0.3186394,0.36013,0.4022501,0.4450965,0.4887764,0.5334097,0.5791322,
    0.626099,0.6744898,0.7245144,0.7764218,0.8305109,0.8871466,0.9467818,
    1.00999,1.077516,1.150349,1.229859,1.318011,1.417797,1.534121,1.67594,
    1.862732,2.153875
};
static double d[31] = {
    0.0,0.0,0.0,0.0,0.0,0.2636843,0.2425085,0.2255674,0.2116342,0.1999243,
    0.1899108,0.1812252,0.1736014,0.1668419,0.1607967,0.1553497,0.1504094,
    0.1459026,0.14177,0.1379632,0.1344418,0.1311722,0.128126,0.1252791,
    0.1226109,0.1201036,0.1177417,0.1155119,0.1134023,0.1114027,0.1095039
};
static double t[31] = {
    7.673828E-4,2.30687E-3,3.860618E-3,5.438454E-3,7.0507E-3,8.708396E-3,
    1.042357E-2,1.220953E-2,1.408125E-2,1.605579E-2,1.81529E-2,2.039573E-2,
    2.281177E-2,2.543407E-2,2.830296E-2,3.146822E-2,3.499233E-2,3.895483E-2,
    4.345878E-2,4.864035E-2,5.468334E-2,6.184222E-2,7.047983E-2,8.113195E-2,
    9.462444E-2,0.1123001,0.136498,0.1716886,0.2276241,0.330498,0.5847031
};
static double h[31] = {
    3.920617E-2,3.932705E-2,3.951E-2,3.975703E-2,4.007093E-2,4.045533E-2,
    4.091481E-2,4.145507E-2,4.208311E-2,4.280748E-2,4.363863E-2,4.458932E-2,
    4.567523E-2,4.691571E-2,4.833487E-2,4.996298E-2,5.183859E-2,5.401138E-2,
    5.654656E-2,5.95313E-2,6.308489E-2,6.737503E-2,7.264544E-2,7.926471E-2,
    8.781922E-2,9.930398E-2,0.11556,0.1404344,0.1836142,0.2790016,0.7010474
};
static long inr;
static double snorm,u,ss,ustar,anr,w,y,tt;
    u = rnd();   /* was ranf--JKP*/
    ss = 0.0;
    if(u > 0.5) ss = 1.0;
    u += (u-ss);
    u = 32.0*u;
    inr = (long) (u);
    if(inr == 32) inr = 31;
    if(inr == 0) goto S100;
/*
                                START CENTER
*/
    ustar = u-(double)inr;
    anr = *(a+inr-1);
S40:
    if(ustar <= *(t+inr-1)) goto S60;
    w = (ustar-*(t+inr-1))**(h+inr-1);
S50:
/*
                                EXIT   (BOTH CASES)
*/
    y = anr+w;
    snorm = y;
    if(ss == 1.0) snorm = -y;
    return snorm;
S60:
/*
                                CENTER CONTINUED
*/
    u = rnd();                /*was ranf--JKP*/
    w = u*(*(a+inr)-anr);
    tt = (0.5*w+anr)*w;
    goto S80;
S70:
    tt = u;
    ustar = rnd();                /*was ranf--JKP*/
S80:
    if(ustar > tt) goto S50;
    u = rnd();               /*was ranf--JKP*/
    if(ustar >= u) goto S70;
    ustar = rnd();               /*was ranf--JKP*/
    goto S40;
S100:
/*
                                START TAIL
*/
    inr = 6;
    anr = *(a+31);
    goto S120;
S110:
    anr += *(d+inr-1);
    inr += 1;
S120:
    u += u;
    if(u < 1.0) goto S110;
    u -= 1.0;
S140:
    w = u**(d+inr-1);
    tt = (0.5*w+anr)*w;
    goto S160;
S150:
    tt = u;
S160:
    ustar = rnd();               /*was ranf--JKP*/
    if(ustar > tt) goto S50;
    u = rnd();               /*was ranf--JKP*/
    if(ustar >= u) goto S150;
    u = rnd();               /*was ranf--JKP*/
    goto S140;
}

/*----generate real number from N(mu, sd) (by STRUCTURE)-------------------------------------------------------------------------*/
double RNormal(double mu, double sd) 
     /* Returns Normal rv with mean mu, variance sigsq.   
        Uses snorm function of Brown and Lovato.  By JKP*/
{
  return (mu + sd*snorm());
}

/*----Normal density--------------------------------------------------------------------------------------------------------------*/
double LogDNormal(double x, double mu, double sd)
{
	return(-0.9189385 - 0.5*log(sd) - 0.5*pow(x-mu,2.0)/pow(sd,2.0));/* -0.5*log(2pi) = -0.9189385 */
}

/*----calculate log of the gamma function (from STRUCTURE)-------------------------------------------------------------------------*/
double mylgamma(double zpo)
{
/* LGAMMA function

   double_value = lgamma(<double_value > 0.>)

   returns the natural log of the gamma function

Uses Lanczos-type approximation to ln(gamma) for zpo > 0.   
Reference:                                          
 Lanczos, C. 'A precision approximation of the gamma   
    function', J. SIAM Numer. Anal., B, 1, 86-96, 1964.   

Original was in FORTRAN
Accuracy: About 14 significant digits except for small regions   
          in the vicinity of 1 and 2.  
Programmer: Alan Miller   
          CSIRO Division of Mathematics & Statistics   
Latest revision - 17 April 1988   
  
Translated and modified into C by Peter Beerli 1997 
Tested against Mathematica's Log[Gamma[x]]
*/

  double alg[9] = { 0.9999999999995183, 676.5203681218835,
    -1259.139216722289, 771.3234287757674, -176.6150291498386,
    12.50734324009056, -0.1385710331296526, 9.934937113930748e-6,
    1.659470187408462e-7 };
  double lnsqrt2pi = 0.9189385332046727;
  double result;
  long jlg;
  double tmp;
  if (zpo <= 0.)
    {
      fprintf(stderr,"lgamma function failed with wrong input (%lf)\n",zpo);
      exit(-1);
    }
  result = 0.;
  tmp = zpo + 7.;
  for (jlg = 9; jlg >= 2; --jlg)
    {
      result += alg[jlg - 1] / tmp;
      tmp -= 1.;
    }
  result += alg[0];
  result = log (result) + lnsqrt2pi - (zpo + 6.5) + (zpo - 0.5) * log (zpo + 6.5);
  return result;
}

/*----Density of the inverse chi square distribution------------------------------------------------------------------------------*/
double LogDInchi(double x, double nu, double s2)
{
	return(0.5*nu*log(0.5*nu)-mylgamma(0.5*nu)+0.5*nu*log(s2)-(0.5*nu+1.0)*log(x)-0.5*nu*s2/x);
}

/*----generate a real number from inverse chi square distribution, chi ( v, s2 ) ------------------------------------------------------*/
double Rinchi ( double v, double s2 )
{
	return ( v * s2 / Rchi ( v ) );
}

/*----LU decomposition of symmetric matrix----------------------------------------------------------------------------------------*/
/* return determinant */
double LUdesym (const double *sigma, double *low, int n)
{
	int	ii, jj, kk, iinjj;
	double det;

	det = 1.0;

	low [0] = sqrt(sigma [0]);

	for (ii=0; ii<n; ii++)
		for (jj=0; jj<n; jj++)
		{
			if (ii>0 || jj>0)
				if (ii>=jj)	low [ii*n + jj] = sigma [ii*n + jj];
		}

	for (ii=1; ii<n; ii++)
		low [ii*n] /= low [0];


	for (jj=1; jj<n; jj++)
		for (ii=jj; ii<n; ii++)
		{
			iinjj  = ii*n + jj;
			if (ii==jj)	/* calculate the diagonal elements */
			{
				for (kk=0; kk<jj; kk++)
					low[iinjj] -= pow(low[ii*n + kk], 2);

				low[iinjj] = sqrt(low[iinjj]);

/**/
						if(low[iinjj]<0.0)
						{
							printf("jj%d %f\n", jj, low[iinjj]);
							for(kk=0; kk<jj; kk++)
								printf("kk%d %f\n", kk, low[ii*n + kk]);

						}
/**/				
			}
			else		/* calculate the off-diagonal elements */
			{
				for (kk=0; kk<jj; kk++)
					low[iinjj] -= (low[ii*n + kk] * low [jj*n + kk]);
				low[iinjj] /=  low[jj*n + jj];
			}

		}



	for (ii=0;ii<n;ii++)
		det *= pow(low[ii*n+ii],2);

	return(det);

}

/*----Multiply square matrices. B = A %*% B %*% t(A)-------------------------------------------------------------------------------*/
int MultiplySquareMatrices (const double *A, double *B, int dim)
{
	int ii,jj,kk, Pos;

	for (ii=0; ii<dim; ii++)
		for (jj=0; jj<dim; jj++ )
		{
			Pos = ii*dim + jj;
			Intermediate [Pos] = 0.0;
			for (kk=0; kk<dim; kk++ )
				Intermediate [Pos] += (A [ii*dim + kk] * B [kk*dim + jj]);
		}

	for (ii=0; ii<dim; ii++)
		for (jj=0; jj<dim; jj++ )
		{
			Pos = ii*dim + jj;
			B [Pos] = 0.0;
			for (kk=0; kk<dim; kk++ )
				B [Pos] += (Intermediate [ii*dim + kk] * A [jj*dim + kk]);
		}
			
	return(0);
}

/*----generate random values from Wishart distribution, W (def, sigma), 
	where def is the degree of freedom and sigma is a symmetric matrix---------------------------------------------------------------*/
int RWishart (const double *sigma, double def, int dim, double *sample)
{
	int ii, jj, kk, mm, Ne;
	double temp;
	
	Ne = dim*dim;
	for (ii=0; ii<Ne; ii++)
		Lowertri[ii] = 0.0;

	temp = LUdesym ( sigma, Lowertri, dim);

/*

			printf ("Lowertri\n");
			for (ii=0; ii<dim; ii++ )
			{
				printf ( "%f ", Lowertri[ii*dim+ii] );
			}
			putchar ('\n');
*/

	for (jj=0; jj<dim; jj++)
	{
		Vi [jj] = Rchi ( def - (double)jj);
		sqrtVi [jj] = sqrt (Vi[jj]);
		for (kk=jj+1; kk<dim; kk++)
				Nij [jj*dim + kk] = snorm();
	}

	sample [0] = Vi[0];
	for (jj=1; jj<dim; jj++)
	{
		sample [dim*jj + jj] = Vi[jj];
		for (kk=0; kk<jj; kk++)
			sample [jj*dim + jj] += pow (Nij [kk*dim + jj], 2);

	}

	for (jj=1; jj<dim; jj++)
	{
		sample [jj] = Nij [jj] * sqrtVi[0];
		sample [jj*dim] = sample [jj]; 
	}

	for (jj=1; jj<dim; jj++)
	{
		for (kk=jj+1; kk<dim; kk++)
		{
			sample [jj*dim + kk] = Nij [jj*dim + kk] * sqrtVi[jj];				
			for (mm=0; mm<jj; mm++)
				sample [jj*dim + kk] += (Nij [mm*dim + jj] * Nij [mm*dim + kk]);

			sample [kk*dim + jj] = sample [jj*dim + kk];
		}
	}

	MultiplySquareMatrices ( Lowertri, sample, dim);

	return (0);
}

/*----Indicator for symmetric matrix---------------------------------------------------------------------------------------------------*/
/* Position in a symmetric matrix, [a, b], is converted to the corresponding position of the array 
that consists of the upper triangular part of the symmetric matrix*/
int	PositionInSym (int a, int b, int dim)
{
	if (b>=a) { return((int) (a*dim - (a-1)*a/2 + b - a ) ); } else { return ( (int) (b*dim - (b-1)*b/2 + a - b) ); }
}

/*----Multiply a vector with a symmetric matrix. result = vector %*% matrix ------------------------------------------------------------*/
int MultiplyVectorSymMatrixLstruct (Lstruct *vector, int tr, double *matrix, double *result, int dim)
{
	int	ii,jj;
	double temp;

	for (ii=0; ii<dim; ii++ )
	{
		temp=0.0;
		for (jj=0; jj<dim; jj++ )
			temp += (vector[tr*dim + jj].currenteffect * matrix [ PositionInSym(jj,ii,dim) ]);

		result [ii] = temp;
	}

	return(0);

}

/*----Multiply vectors. result = A %*% B ------------------------------------------------------------------------------------*/
int MultiplyVectors (const double *A, const double *B, double *result, int dim)
{
	int ii;

	result [0] = 0.0;
	for (ii=0; ii<dim; ii++)
		result [0] += A[ii]*B[ii];

	return(0);
}

/*----Multiply vectors. result = A %*% B ------------------------------------------------------------------------------------*/
int MultiplyVectorsLstruct (const double *A, const Lstruct *B, int tr, double *result, int dim)
{
	int ii;

	result [0] = 0.0;
	for (ii=0; ii<dim; ii++)
		result [0] += A[ii]*B[tr*dim + ii].currenteffect;

	return(0);
}

/*----perform Gauss-Jordan algorithm 
	from GaussJordan.exe source cord-----------------------------------------------------------------------------------------*/
int GaussJordan (double *A, double *result, int dim)
{
	int		ii,jj, kk, Ne;
	int		Pivotrow, Currentrow, Treatingrow;
	double	Max, temp, invDiag;

	Ne = dim*dim;

	for (ii=0; ii<dim; ii++)
		for (jj=0; jj<dim; jj++)
			if ( ii==jj ) { result [ii*dim + jj] = 1.0; } else { result [ii*dim + jj] = 0.0;}

	for (ii=0; ii<dim; ii++)
	{
		Max = 0.0;
		for (jj=ii; jj<dim; jj++)
		{
			temp = fabs(A[jj*dim + ii]);
			if (temp > Max)
			{
				Max = temp;
				Pivotrow = jj;
			}
		}
/**/
		if (Max == 0.0)
		{
			printf ( "Inversion error\n" );
/*			for(jj=ii; jj<dim; jj++)
				printf ("jj%d %f\n", jj+1, A[jj*dim + ii]);*/
		}
/**/
		Currentrow = ii*dim;

		if ( Pivotrow != ii )
		{
			Pivotrow *= dim;
			for (jj=0; jj<dim; jj++)
			{
				temp = A [Currentrow + jj];
				A [Currentrow + jj] = A [Pivotrow + jj];
				A [Pivotrow + jj] = temp;

				temp = result [Currentrow + jj];
				result [Currentrow + jj] = result [Pivotrow + jj];
				result [Pivotrow + jj] = temp;
			}
		}

		invDiag = 1.0/A[Currentrow + ii];

		for (jj=0; jj<dim; jj++)
		{
			A[Currentrow + jj] *= invDiag;
			result[Currentrow + jj] *= invDiag;
		}

		for (jj=0; jj<dim; jj++)
		{
			if (jj!=ii)
			{
				Treatingrow = jj*dim;
				temp = A[Treatingrow + ii];

				for (kk=0; kk<dim; kk++)
				{
					A[Treatingrow + kk] -= (temp * A[Currentrow + kk]);
					result[Treatingrow + kk] -= (temp * result[Currentrow + kk]);
/*
					if(jj==kk)
						if (result[Treatingrow + kk] < 1.0e-20)
							result[Treatingrow + kk] = 1.0e-20;
*/
				}
			}
		}
	}		/* ii */

	return(0);
}


/*----Returns a random number between 0 and total-1, according the probalities list, *probs.-------------------------------------*/ 
    /*sum is the sum of the probabirities (by STRUCTURE with some modifications) */
int PickAnOption(int total, double sum, const double *probs)
{
  int option;
  double random;
  double sumsofar = 0.0;

  random = RandomReal(0,sum);     /*Get uniform random real in this range*/
  for (option=0; option<total; option++) /*Figure out which popn this is*/
    {
      sumsofar += probs[option];
	if (random <= sumsofar) break;
    }

  return option;

}

/*----Returns a random number between 0 and total-1, according the probalities list, *probs.-------------------------------------*/ 
	/* probs is log */
int PickAnOptionL(int total, double *probs)
{
  int option;
  double random, sum, max, value;
  double sumsofar = 0.0;

  max = probs [0];
  if ( max = 0.0 ) max = probs [1];
  for (option=1; option<total; option++ )
  {
	  value = probs [option];
	  if ( value > max && value != 0.0 )
		  max = value;
  }

  sum = 0.0;
  for (option=0; option<total; option++ )
  {
	  value = probs [option];
	  if ( value != 0.0 ) 
	  {
		  value -= max;
		  probs [option] = exp ( value );
	  }
	  sum += probs [option];
  }

  random = RandomReal(0,sum);     /*Get uniform random real in this range*/
  for (option=0; option<total; option++) /*Figure out which popn this is*/
    {
      sumsofar += probs[option];
	if (random <= sumsofar) break;
    }

  return option;

}

/*----Initialize classes according to the Dirichlet process prior ----------------------------------------------------------------*/
int Initializeclass (int *lambdaclass, double alpha, int N)
{
	int ii, Nc;
	double *Nic;
	
	Nic = (double*) calloc (N, sizeof(double) );
	Nic [0] = 1.0;
	Nc = 1;
	
	for (ii=1; ii<N; ii++)
	{
		Nic [Nc] = alpha;
		lambdaclass[ii] = PickAnOption (Nc+1, (double)ii+alpha, Nic);

		if (lambdaclass[ii] == Nc){
			Nic [Nc] = 1.0;
			Nc ++;
		}
		else
		{
			Nic[Nc] -= alpha;
			Nic [lambdaclass[ii]] ++;
		}

	}

	free (Nic);

	return(Nc);

}

/* return the exponential term of the likelihood function -------------------------------------------------------------------*/
double LikelihoodExp (int K, double *iR, double *Eri, double *Products )
{
	int ii, jj;
	double temp;

	for (ii=0; ii<K; ii++)
	{
		Products [ii] = 0.0;
		for (jj=0; jj<K; jj++)
			Products [ii] += ( Eri [jj] * iR [ jj*K + ii]);
	}

	temp = 0.0;
	for (ii=0; ii<K; ii++)
		temp += (Products [ii] *  Eri [ii] );

	temp *= (-0.5);

	return(temp);
}

/*----alert when allocation of memory is failed-----------------------------------------------------------------------------------*/
void alert ( void )
{
	puts ( "memory allocation failure. try again" );
}

/* de Boor Cox algorithm for B spline----------------------------------------------------------------------------------------------*/
double deBoorCox (double t, int i, double t0, double a, int n)
{
	/* t: data point 
	   i: ith control point (0<=i<Nsp)
	   t0:first knot
	   a: interval of knots
	   n: order of spline
	*/
	double ti, ti1, result=0.0;

	ti=t0 + a*(double)i;
	ti1=ti + a;
	if(n==0){
		if(ti<=t&&t<ti1) result=1.0; else result=0.0;
	}else{
		result=(t-ti)/((double)n*a)*deBoorCox(t,i,t0,a,n-1)+(ti+a*(double)(n+1)-t)/((double)n*a)*deBoorCox(t,i+1,t0,a,n-1);
	}
	return(result);
}

double deBoorCox2 (double x, int i, double *knot, int m)
{
	/* t: data point 
	   i: ith spline (0<=i<Nsp)
	   knot: positions of knots
	   m: order of spline (1<=m<=M)
	*/
	double ti, ti1, result=0.0, temp1, temp2;

	ti=knot[i];
	ti1=knot[i+1];
	if(m==1){
		if(ti<=x&&x<ti1) result=1.0; else result=0.0;
	}else{
		if((knot[i+m-1]-ti)>0.0) temp1 = (x-ti)/(knot[i+m-1]-ti)*deBoorCox2(x,i,knot,m-1); else temp1 = 0.0;
		if((knot[i+m]-ti1)>0.0) temp2 = (knot[i+m]-x)/(knot[i+m]-ti1)*deBoorCox2(x,i+1,knot,m-1); else temp2 = 0.0;
		result= temp1 + temp2;
	}
	return(result);
}


/*----Multiple trait analysis------------------------------------------------------------------------------------------------------*/
int mta05 (double Missingvalue, int K, int P, int N, double va, double ve, int Burnin, int Ite, int Thi, int CovR, int *Ft, int F, int *cumFt,
			Ystruct *Y, Lambdastruct *L, Lstruct *X, Lstruct *Z, int SEM, int M, double Lambda0, double Lambdasigmafixed, 
			int	Polygenic, double *iA, double *Va, double *Ve, double *ReCor, double *Nu, double *S2, int Identitymatrix, double *Mean, double *Sd,
			double *SamplediG, double *SamplediR, double *SampledLikelihood, int SplinePrior, int AddIntercept)
{
	/*
	===CAUTION======================
		The likelihood function used here is
			
			(2pi)^(-KN/2) * det(R)^(-N/2) * exp ( - 0.5* (y'-XB-Zu)' iR (y'-XB-Zu) )

		where y' indicates y adjusted for the recursive effects. This function is valid when only the recursive effects are included in the model.
		If simultaneous effects are included, this function should be multiplied by det(Lambda).
	==================================
	*/

	/* MCMC iteration */
	int		Totalite, mcmc;

	/* For repeat statement */
	int		ii, jj, mm, nn;
	int		trait, record, additive, fixed, cla;

	/* Generic use */
	double	temp, temp2, prop;
	int		target;

	/* Assignment of individuals. Used only when the Diriclet process is applied */
	int		Naux = 4;		/* number of auxilialy parameters is set to 4 */		
	int		Arraysize, Originalclass, UpdateEr;
	double	*Auxvalues, *Probs, *Copydouble, Originallambda;
	int		*Copyint;

	/* posterior variance */
	double	iC;

	/* residual (co)variances */
	double	*R;

	/*inverted genetic (co)variances, inverted environmental (co)variances */
	double	*iG, *iR;

	/* Sa, Se, inverted Sa, inverted Se*/
	double	*Sa, *Se, *iSa, *iSe;

	/* For calculation of Sa */
	double	*Saproducts;

	/* For Metropolis update */
	double	MHprob, *propR, *propiR, DetpropR, propPm;
	double	SdRprop=0.2;		/* [SD of observations] * SdRprop is used for SD of proposal normal distribution of resiudal variance (when CovR=1) */
	double	SdPm=0.2;
	double	BecomeVague=1000; /* Prior variances of P1 and P2 are (BecomeVague * L[mm].sigmaPi) */
	double	sumPm2;				/* sum of P^2 used to update sigmaPi */
	int		NacceptedR=0;
	int		UseVaguePrior;		/* The first "UseVaguePrior" basis splines use a vague prior. Usually 1, but 2 if the intercept is included. */

	/* For sampling */
		/* Indicator */
		int		DoSampling, AddtoMean;

		/* For sampling parameters */
		int		Ns, CumNs;

	/* For log likelihood calculation */
	double	*Dummy, *Products, *Eri, DetR, LikeConstant;

	/* Products and sums of hyperparameters */
	double  veN, vaP;

	/* Products of varables */
	int		KN, KK, PP12, KK12, PK, NcNaux, Naux1, NcNaux1, effectpos, diagonal;
	double	P2;

	/* Products of parameters */
	KN = K*N;
	KK = K*K;
	PK = P*K;
	PP12 = P*(P+1)/2;
	KK12 = K*(K+1)/2;
	Naux1 = Naux - 1;

	veN = ve + (double)N;
	vaP = va + (double)P;

	P2 = (double) P - 2.0;	/* to generate random values from the prosterior distribution of genetic variances. Used when polygenic==1 and K==1 */

	/* Total number of MCMC iterations */
	Totalite = Burnin + Ite;

    /* Number of samples */
    Ns = Totalite / Thi;

	/* For likelihood calculation */
	Dummy = (double*) calloc ( KK, sizeof(double));								if(Dummy==NULL) alert();
	Products = (double*) calloc ( K, sizeof(double));							if(Products==NULL) alert();
	Eri = (double*) calloc ( K, sizeof(double));								if(Eri==NULL) alert();
	LikeConstant = pow(2.0 * Pi,-0.5*(double)K);

	/* For SEM, records are made missing when their upstream or downstream traits are missed */
	for (mm=0; mm<M; mm++)
		for (record=0; record<N; record++)
		{
			if(Y[L[mm].up].observations[record] == Missingvalue)
				Y[L[mm].down].observations[record] = Missingvalue;
			if(Y[L[mm].down].observations[record] == Missingvalue)
				Y[L[mm].up].observations[record] = Missingvalue;
		}

	/* Check missing values */
	for (trait=0; trait<K; trait++)
	{
		for ( record=0; record<N; record++ )
			if (Y[trait].observations[record]==Missingvalue)
				Y[trait].nmiss ++;

		if(Y[trait].nmiss > 0)
		{
			Y[trait].whichrec = (int*) calloc (Y[trait].nmiss, sizeof(int));	if(Y[trait].whichrec==NULL) alert();
			nn = 0;
			for (record=0; record<N; record++)
				if (Y[trait].observations[record]==Missingvalue) { Y[trait].whichrec[nn] = record; nn++;}
		}
	}

	/* Mean and SD of observations */
	for (trait=0, Mean[trait]=0.0, Sd[trait]=0.0; trait<K; trait++)
	{
		for(record=0; record<N; record++)
			if(Y[trait].observations[record] != Missingvalue)
				Mean[trait] += Y[trait].observations[record];

		Mean[trait] /= (double) (N - Y[trait].nmiss);

		for(record=0; record<N; record++)
			if(Y[trait].observations[record] != Missingvalue)
				Sd[trait] += pow( (Y[trait].observations[record] - Mean[trait]), 2.0);

		Sd[trait] /= (double) (N - Y[trait].nmiss - 1);
		Sd[trait] = sqrt(Sd[trait]);

		printf ("trait %d Mean %f Sd %f\n", trait, Mean[trait], Sd[trait]);

		Y[trait].stobs = (double*) calloc (N, sizeof(double));
		memcpy (Y[trait].stobs, Y[trait].observations, sizeof(double)*N);

		/* Quit centering 
		for (record=0; record<N; record++)
			if(Y[trait].observations[record] != Missingvalue)
				Y[trait].stobs [record] = (Y[trait].observations [record] - Mean[trait]);//Sd[trait]; Quit dividing by SD.
		*/
		/* From here, "stobs" is used for calculation. "observations" is used for missing value recognition */
	}

	/* Objects related with X */
	if(F>0)
	{
		for (trait=0; trait<K; trait++)
		{
			for (fixed=0; fixed<Ft[trait]; fixed++)
			{
				X[cumFt[trait] + fixed].nrec = 0;
				X[cumFt[trait] + fixed].x2 = 0.0;

				for (record=0; record<N; record++)
					if ( X[cumFt[trait] + fixed].covariates[record] != 0.0)
						X[cumFt[trait] + fixed].nrec ++;

				if(X[cumFt[trait] + fixed].nrec > 0)
				{
					X[cumFt[trait] + fixed].whichrec = (int*) calloc (X[cumFt[trait] + fixed].nrec, sizeof (int));	if (X[fixed].whichrec == NULL) alert();

					ii=0;
					for (record=0; record<N; record++)
						if ( X[cumFt[trait] + fixed].covariates[record] != 0.0)
						{
							X[cumFt[trait] + fixed].whichrec [ii] = record;
							X[cumFt[trait] + fixed].x2 += pow (X[cumFt[trait] + fixed].covariates[record], 2.0);
							ii++;
						}
				}
				else
				{
					X[cumFt[trait] + fixed].whichrec = (int*) calloc (1, sizeof(int));
					X[cumFt[trait] + fixed].whichrec[0] = -1;
				}
			}
		}
	}

	/* Objects related with Z */
	if (Polygenic)
	{
		if (Identitymatrix) /* identity matrix. P=N */
		{
			for (trait=0; trait<K; trait++)
				for (additive=0; additive<P; additive++ )
				{
					Z[trait*P + additive].nrec = 1;
					Z[trait*P + additive].x2 = 1.0;
					Z[trait*P + additive].whichrec = (int*) calloc (1, sizeof(int));
					Z[trait*P + additive].whichrec[0] = additive;
				}
		}
		else
		{
			for (trait=0; trait<K; trait++)
			{
				for (additive=0; additive<P; additive++)
				{
					Z[trait*P + additive].nrec = 0;
					Z[trait*P + additive].x2 = 0.0;

					for (record=0; record<N; record++)
						if (Z[additive].covariates[record]>0.0 )
							Z[trait*P + additive].nrec ++;

					if (Z[trait*P + additive].nrec > 0)
					{
						Z[trait*P + additive].whichrec = (int*) calloc (Z[trait*P + additive].nrec, sizeof(int));
						nn=0;
						for (record=0; record<N; record++)
						{
							if (Z[additive].covariates[record]>0.0)
							{
								Z[trait*P + additive].whichrec [nn] = record;
								Z[trait*P + additive].x2 += pow (Z[additive].covariates[record], 2.0);
								nn ++;
							}
						}
					}
					else
					{
						Z[trait*P + additive].whichrec = (int*) calloc (1, sizeof(int));
						Z[trait*P + additive].whichrec [0] = -1; 
					}
				}
			}
		}		/* if (Zfile[0] == '0') */
	}			/* Polygenic*/

	/* Initialization */
	if (Polygenic)
	{
		iG = (double*) calloc ( KK, sizeof (double) );
		for (trait=0; trait<K; trait++)
			iG[trait*K + trait] = 0.5 * Sd[trait] *Sd[trait];
		Sa  = (double*) calloc (KK, sizeof (double));
		iSa = (double*) calloc (KK, sizeof (double));
		Saproducts = (double*) calloc ( P, sizeof (double) );						if(Saproducts==NULL) alert();

		for (trait=0; trait<K; trait++)
		{
			temp = sqrt( 1.0/iG [trait*K + trait]);
			for (additive=0; additive<P; additive++)
				Z[trait*P + additive].currenteffect = RNormal ( 0.0, temp);
		}
	}

	R = (double*) calloc ( KK, sizeof (double) );		if(R==NULL) alert();
	iR =  (double*) calloc ( KK, sizeof (double) );		if(iR==NULL) alert();
	if(CovR==1)
	{
		propR  =  (double*) calloc ( KK, sizeof (double) );		if( propR==NULL) alert();
		propiR  =  (double*) calloc ( KK, sizeof (double) );	if( propiR==NULL) alert();
		for (trait=0; trait<K; trait++)
		{
			for(ii=trait; ii<K; ii++)
			{
				R [trait*K + ii] = 0.5 * ReCor [trait*K + ii] * Sd[trait] * Sd[ii];
				R [ii*K + trait] = R [trait*K + ii];
			}
		}
		memcpy(Dummy,R,sizeof(double)*KK);
		DetR = LUdesym ( Dummy, iR, K);
		GaussJordan ( Dummy, iR, K);
	}
	else
	{
		for (trait=0, DetR=1.0; trait<K; trait++)
		{
			nn = trait*K + trait;
			R [nn] = 0.5 * Sd[trait] * Sd[trait];
			iR [nn] = 1.0/R [nn];
			DetR *= R [nn];
		}
	}
	Se =  (double*) calloc ( KK, sizeof (double) );		if(Se==NULL) alert();
	iSe = (double*) calloc ( KK, sizeof (double) );		if(iSe==NULL) alert();

	if (SEM) 
	{
		for (mm=0; mm<M; mm++)
		{
			switch (L[mm].update) {

			case 1:
				L[mm].lambdasigma = Lambdasigmafixed;
				L[mm].meanlambda = Lambda0;
				L[mm].lambdahat = (double*) calloc (1, sizeof(double));
				L[mm].vlm = (double*) calloc (1, sizeof(double));

				L[mm].y2 = (double*) calloc (1, sizeof(double));
				for (record=0; record<N; record++)
					if (Y[L[mm].up].observations [record] != Missingvalue)
						L[mm].y2[0] += pow (Y[L[mm].up].stobs [record], 2);

				L[mm].lambda = (double*) calloc ( N, sizeof(double));
				temp = RandomReal (-1.0, 1.0);										/* the initial values of Lambda are generated from Uniform (-1.0, 1.0) */
				for (record=0; record<N; record++) L[mm].lambda[record] = temp;
				break;

			case 2:
				L[mm].lambdaclass = (int*) calloc (N, sizeof(int));
				L[mm].nc = Initializeclass ( L[mm].lambdaclass, L[mm].alpha, N);

				L[mm].nic = (int*) calloc (L[mm].nc, sizeof(int));
				for (record=0; record<N; record++)
					if (Y[L[mm].up].observations [record] != Missingvalue)
						L[mm].nic [L[mm].lambdaclass[record]] ++;

				L[mm].lambda = (double*) calloc ( N, sizeof(double));
				temp = RandomReal (-1.0, 1.0);										/* the initial values of Lambda are generated from Uniform (-1.0, 1.0) */
				for (record=0; record<N; record++) L[mm].lambda[record] = temp;		/* initial value is shared among replicates even if the heterogeneity is assumed */

				L[mm].lambdalist = (double*) calloc (L[mm].nc, sizeof(double));
				for (cla=0; cla<L[mm].nc; cla++)
					L[mm].lambdalist[cla] = temp;

				L[mm].y2 = (double*) calloc (N, sizeof(double));
				for (record=0; record<N; record++)
					if (Y[L[mm].up].observations [record] != Missingvalue)
						L[mm].y2[record] = pow( Y [ L[mm].up].stobs[record], 2);

				Auxvalues = (double*) calloc ( Naux, sizeof(double) );
				break;

			case 3:
				/*search maximum and minimum observed values */
				temp=1e-100; temp2=1e+100;
				for (record=0; record<N; record++)
				{
					if(Y[L[mm].up].observations [record] != Missingvalue)
					{
						if(Y[L[mm].up].stobs[record]>temp) temp=Y[L[mm].up].stobs[record];
						if(Y[L[mm].up].stobs[record]<temp2) temp2=Y[L[mm].up].stobs[record];
					}
				}

				/* ===========================================================================
				Scripts when the first date point is located at the duplicated knots
				L[mm].a = (temp-temp2)/(double)(L[mm].Nsp-L[mm].Ord+1);
				L[mm].knot = (double*)calloc(L[mm].Nsp+L[mm].Ord, sizeof(double));
				for(ii=0; ii<L[mm].Ord; ii++)
				{
					L[mm].knot[ii] = temp2; L[mm].knot[L[mm].Nsp+L[mm].Ord-1-ii] = temp;
				}
				for(ii=L[mm].Ord; ii<L[mm].Nsp; ii++)
				{
					L[mm].knot[ii] = L[mm].knot[ii-1] + L[mm].a;
				}
				==========================================================================*/

				/* ===========================================================================
				Scripts when the first date point is located between the duplicated knots and the next knot*/
				L[mm].a = (temp-temp2)/(double)(L[mm].Nsp-L[mm].Ord);
				L[mm].knot = (double*)calloc(L[mm].Nsp+L[mm].Ord, sizeof(double));
				for(ii=0; ii<L[mm].Ord; ii++)
				{
					L[mm].knot[ii] = temp2-0.5*L[mm].a; L[mm].knot[L[mm].Nsp+L[mm].Ord-1-ii] = temp+0.5*L[mm].a;
				}
				for(ii=L[mm].Ord; ii<L[mm].Nsp; ii++)
				{
					L[mm].knot[ii] = L[mm].knot[ii-1] + L[mm].a;
				}
				/*==========================================================================*/

				/* ===========================================================================
				Scripts when the first date point is located at the next knot of the duplicated knots
				L[mm].a = (temp-temp2)/(double)(L[mm].Nsp-L[mm].Ord-1);
				L[mm].knot = (double*)calloc(L[mm].Nsp+L[mm].Ord, sizeof(double));
				for(ii=0; ii<L[mm].Ord; ii++)
				{
					L[mm].knot[ii] = temp2-L[mm].a ; L[mm].knot[L[mm].Nsp+L[mm].Ord-1-ii] = temp+L[mm].a;
				}
				for(ii=L[mm].Ord; ii<L[mm].Nsp; ii++)
				{
					L[mm].knot[ii] = L[mm].knot[ii-1] + L[mm].a;
				}
				===============================================================================*/

				/*---
				for(ii=0; ii<(L[mm].Nsp+L[mm].Ord); ii++)
				{
					printf("%f ",L[mm].knot[ii]);
				}
				printf("\n");
				---*/

				/* Prepare objects related with B */
				L[mm].B = (double*) calloc ( L[mm].Nsp*N, sizeof(double));
				L[mm].HaveithSP = (int*) calloc ( N*L[mm].Nsp, sizeof(int));
				L[mm].NhaveithSP = (int*) calloc ( L[mm].Nsp, sizeof(int));
				for(record=0; record<N; record++)
					if(Y[L[mm].up].observations [record] != Missingvalue)
						for(ii=0; ii<L[mm].Nsp; ii++)
							L[mm].B[record*L[mm].Nsp+ii] = deBoorCox2(Y[L[mm].up].stobs[record],ii,L[mm].knot,L[mm].Ord);
				if(AddIntercept) /* the first basis spline is replaced by the intercept */
					for(record=0; record<N; record++)
						if(Y[L[mm].up].observations [record] != Missingvalue)
							L[mm].B[record*L[mm].Nsp+0] = 1.0;				

				/*---
				record=0;
				for(ii=0; ii<L[mm].Nsp; ii++)
				{
					printf("%f\n", L[mm].B[record*L[mm].Nsp+ii]);
				}
				record=4;
				for(ii=0; ii<L[mm].Nsp; ii++)
				{
					printf("%f\n", L[mm].B[record*L[mm].Nsp+ii]);
				}
				---*/

				for(ii=0; ii<L[mm].Nsp; ii++)
					for(record=0; record<N; record++)
						if((Y[L[mm].up].observations [record] != Missingvalue)&&(L[mm].B[record*L[mm].Nsp+ii]>0.0))
						{
							L[mm].HaveithSP[ii*N+L[mm].NhaveithSP[ii]] = record;
							L[mm].NhaveithSP[ii]++;
						}

				/* Initialize P. The initial values are generated from Uniform (-1.0, 1.0) */
				L[mm].P = (double*) calloc (L[mm].Nsp, sizeof(double));
				for(ii=0; ii<L[mm].Nsp; ii++) L[mm].P[ii] = RandomReal (-1.0, 1.0);
				L[mm].lambda = (double*) calloc ( N, sizeof(double));
				for (record=0; record<N; record++) 
				{
					for(ii=0; ii<L[mm].Nsp; ii++)
						L[mm].lambda[record] += L[mm].P[ii] * L[mm].B[record*L[mm].Nsp+ii];
				}
				/* Initialize sigmaPi */
				L[mm].sigmaPi = 1.0;
				/* Define prior distributions */
				if(SplinePrior==1)
				{
					if(AddIntercept) UseVaguePrior=2; else UseVaguePrior=1;
				}
				else
				{
					if(AddIntercept) UseVaguePrior=3; else UseVaguePrior=2;
				}
				break;

			} /* switch (L[mm].update) */
		}	  /* for (mm=0; mm<M; mm++) */
	}

	if (F>0)
		for(trait=0; trait<K; trait++)
			for(fixed=0; fixed<Ft[trait]; fixed++)
				if(fixed==0) X[cumFt[trait]+fixed].currenteffect = Mean[trait]; else X[cumFt[trait]+fixed].currenteffect = 0.0;

	/* Calculate residual errors */
	for (trait=0; trait<K; trait++)
	{
		Y[trait].errors = (double*) calloc (N, sizeof(double));		if(Y[trait].errors==NULL) alert();
		memcpy (Y[trait].errors, Y[trait].stobs, sizeof(double)*N);

		for (ii=0; ii<Y[trait].nmiss; ii++)
			Y[trait].errors [ Y[trait].whichrec[ii] ] = RNormal (0.0, sqrt(1.0/iR[trait*K + trait])); /* errors of missing records are randomly generated*/

		if (F>0)
		{
			for(fixed=0; fixed<Ft[trait]; fixed++)
				for (ii=0; ii<X[cumFt[trait]+fixed].nrec; ii++)
					Y[trait].errors[ X[cumFt[trait]+fixed].whichrec[ii] ] -= 
					(X[cumFt[trait]+fixed].currenteffect * X[cumFt[trait]+fixed].covariates [X[cumFt[trait]+fixed].whichrec[ii]]);
		}

		if (Polygenic)
		{
			for(additive=0; additive<P; additive++)
				for (ii=0; ii<Z[trait*P + additive].nrec; ii++)
					Y[trait].errors[ Z[trait*P + additive].whichrec[ii] ] -=
					(Z[trait*P + additive].currenteffect * Z[additive].covariates[Z[trait*P + additive].whichrec[ii]]);
		}
	}

	if (SEM)
		for(mm=0; mm<M; mm++)
			for (record=0; record<N; record++)
				if (Y[L[mm].up].observations [record] != Missingvalue)
					Y[L[mm].down].errors[record]  -= ( Y [ L[mm].up].stobs [record] * L[mm].lambda [record] );

	CumNs = -1;

	printf ("start MCMC iterations ...\n");
	putchar ('\n');

	/* MCMC iterations */
	for (mcmc=1; mcmc<=Totalite; mcmc++ )
	{
		/* Print mcmc */
		if ( mcmc % 100 == 0 ){printf ( "%d\n", mcmc);}

		/* Sample or Not */
		if ( mcmc % Thi == 0) 
		{
			DoSampling=1; CumNs++;
			if ( mcmc > Burnin) { AddtoMean=1;} else { AddtoMean=0;}
		} 
		else
		{ 
			DoSampling=0; AddtoMean=0;
		}

		/* Update B */	
		if (F>0)
		{						
			for (trait=0; trait<K; trait++ )
			{
				for (fixed=0; fixed<Ft[trait]; fixed++)
				{
					temp = 0.0; effectpos = cumFt[trait] + fixed;
					for (jj = 0; jj<K; jj++)
						for (record=0; record<X[effectpos].nrec; record++)
							temp += iR[trait*K + jj] 
									* X[effectpos].covariates[ X[effectpos].whichrec[record] ] 
									* Y[jj].errors[ X[effectpos].whichrec[record] ];

					iC = 1.0 / (iR[trait*K + trait] * X[effectpos].x2);
					prop = RNormal( iC * temp + X[effectpos].currenteffect,  sqrt(iC) );

					for (record=0; record<X[effectpos].nrec; record++)
						Y[trait].errors[ X[effectpos].whichrec[record] ] += X[effectpos].covariates[X[effectpos].whichrec[record]] *( X[effectpos].currenteffect - prop);

					X[effectpos].currenteffect = prop;

					if (DoSampling)
						X[effectpos].sampledeffect[CumNs] = prop;
				} 
			}
		}

		/* Update U */
		if (Polygenic)
		{
			for (trait=0; trait<K; trait++ )
			{
				diagonal = trait*K + trait;
				for (additive=0; additive<P; additive++ )
				{
					temp = 0.0; effectpos = trait*P + additive;
					for (ii=0; ii<K; ii++ )
					{
						for (record=0; record<Z[effectpos].nrec; record++)
							temp += iR[trait*K + ii] 
									* Z[additive].covariates[Z[effectpos].whichrec[record]] 
									* Y[ii].errors[ Z[effectpos].whichrec[record] ];

						temp2 = 0.0;
						for (jj=0; jj<P; jj++)
							temp2 += Z[ii*P + jj].currenteffect * iA[ PositionInSym(additive, jj, P) ];
						temp2 *= iG [trait*K + ii];

						temp -= temp2;
					}

					temp += iR[diagonal] * Z[effectpos].x2 * Z[effectpos].currenteffect;
					temp += iG[diagonal] * iA[ PositionInSym(additive, additive, P) ] * Z[effectpos].currenteffect;

					iC = 1.0/( iR [diagonal] * Z[effectpos].x2 + iG [diagonal] * iA [ PositionInSym(additive, additive, P) ] );

					prop = RNormal( iC * temp, sqrt(iC) );

					for (record=0; record<Z[effectpos].nrec; record++)
							Y[trait].errors[ Z[effectpos].whichrec[record] ] += Z[additive].covariates[ Z[effectpos].whichrec[record] ]*( Z[effectpos].currenteffect - prop);
						
					Z[effectpos].currenteffect = prop;

					if (DoSampling)
						Z[effectpos].sampledeffect[CumNs] = prop;
				}
			}
		}

		/* Update iR */
		switch(CovR){
		case 0: /* independent residuals */
			for (trait=0, DetR=1.0; trait<K; trait++)
			{
				nn = trait*K + trait;
				MultiplyVectors ( Y[trait].errors, Y[trait].errors, &temp, N );
				R [nn] = Rinchi ( (double)N+Nu[trait], (temp+Nu[trait]*S2[trait])/((double)N+Nu[trait]) ); 
				/* => When missing observations are included, N should be replaced by the number of non-missing values */
				iR [nn] = 1.0 / R [nn];
				DetR *= R[nn];
			}
			break;
		case 1: /* correlation is fixed */
			for (trait=0; trait<K; trait++)
			{
				nn = trait*K + trait;
				memcpy(propR, R, sizeof(double)*KK);
				propR[nn] = RNormal (R[nn],SdRprop*Sd[trait]);
				if(propR[nn]>0)
				{
					for(ii=0; ii<K; ii++)
					{
						if(ii!=trait)
						{
							propR[trait*K+ii] *= sqrt(propR[nn]/R[nn]);
							propR[ii*K+trait] *= sqrt(propR[nn]/R[nn]);
						}
					}
					memcpy(Dummy, propR, sizeof(double)*KK);
					GaussJordan(Dummy, propiR, K);
					DetpropR = LUdesym ( propR, Dummy, K);

					MHprob = 0.0;
					MHprob -= (-0.5 * (double) N * log(DetR));
					MHprob += (-0.5 * (double) N * log(DetpropR));
					for (record=0; record<N; record++)
					{
						for (ii=0; ii<K; ii++)
						{
							Eri[ii] =  Y[ii].errors[record];
						}
						MHprob -= LikelihoodExp(K, iR, Eri, Products);
						MHprob += LikelihoodExp(K, propiR, Eri, Products);
					}

					if(Nu[trait]>0.0&&S2[trait]>0.0) /* when informative priors are assigned */
					{
						MHprob += LogDInchi (propR[nn],Nu[trait],S2[trait]);
						MHprob -= LogDInchi (R[nn],Nu[trait],S2[trait]);
					}

					if(log(rnd())<MHprob)
					{
						memcpy(R,propR,sizeof(double)*KK);
						memcpy(iR,propiR,sizeof(double)*KK);
						DetR = DetpropR;
						NacceptedR++;
					}
				}
			}/* trait*/
			break;
		case 2: /* covariance is inferred */
			for (trait=0; trait<K; trait++ )
			{
				for (ii=trait; ii<K; ii++ )
				{
					nn = trait*K+ii;
					MultiplyVectors ( Y[trait].errors, Y[ii].errors, Se+nn, N );
					Se [nn] += Ve [nn];
					Se [ii*K + trait] = Se [nn];
				}
			}				
			GaussJordan ( Se, iSe, K);
			RWishart ( iSe, veN, K, iR);
			memcpy(Dummy, iR, sizeof(double)*KK);
			GaussJordan (Dummy, R, K);
			DetR = LUdesym ( R, Dummy, K);
			break;
		}

		/* Update iG */
		if (Polygenic)
		{
			for (trait=0; trait<K; trait++ )
			{
				MultiplyVectorSymMatrixLstruct ( Z, trait, iA, Saproducts, P);	
											
				for (ii=trait; ii<K; ii++ )
				{
					nn = trait*K+ii;
					MultiplyVectorsLstruct ( Saproducts, Z, ii, Sa+nn, P);
					Sa[nn] += Va [nn];
					Sa[ii*K + trait] = Sa[nn];
				}
			}
			if ( K==1 )
			{
				iG [0] = 1.0 / Rinchi ( P2, Sa[0]/P2 );
			}
			else
			{
				GaussJordan ( Sa, iSa, K );
				RWishart ( iSa, vaP, K, iG);
			}
		}

		/* Update individual classes (only model 2)*/
		if (SEM)
		{
			for (mm=0; mm<M; mm++)
			{
				switch (L[mm].update) {
				case 2:
					for (record=0; record<N; record++)
					{
						if (Y[L[mm].up].observations [record] != Missingvalue)
						{
							UpdateEr = 0;
							L[mm].nic [L[mm].lambdaclass[record]] --;

							if ( L[mm].nic [L[mm].lambdaclass[record]] == 0 )	/* in the case where record is singleton. */
							{

								NcNaux1 = L[mm].nc + Naux1;
								Probs = (double*) calloc ( NcNaux1, sizeof (double));

								for (cla=0; cla<Naux1; cla++)
									Auxvalues [cla] = RNormal ( L[mm].meanlambda, sqrt ( L[mm].lambdasigma ) );

								for (trait=0; trait<K; trait++)
									Eri[trait] =  Y[trait].errors[record];

								for (cla=0; cla<L[mm].nc; cla++)
								{
									Eri [L[mm].down] += (Y [L[mm].up].stobs [record] * ( L[mm].lambda[record] - L[mm].lambdalist[cla] ));
									Probs [cla] += ( LikelihoodExp (K, iR, Eri, Products) );
									Eri [L[mm].down] = Y[L[mm].down].errors[record];
								}

								for (cla=0; cla<Naux1; cla++)
								{
									Eri [L[mm].down] += (Y [L[mm].up].stobs [record] * ( L[mm].lambda[record] - Auxvalues [cla] ));
									Probs [cla+L[mm].nc] += ( LikelihoodExp (K, iR, Eri, Products) );
									Eri [L[mm].down] = Y[L[mm].down].errors[record];
								}

								temp = log ( L[mm].alpha/(double) Naux );
								for (cla=0; cla<L[mm].nc; cla++)
									if (cla == L[mm].lambdaclass[record]) { Probs [cla] += temp; } else { Probs [cla] += log ( (double) L[mm].nic [cla] );}

								for (cla=0; cla<Naux1; cla++)
									Probs [cla+L[mm].nc] += temp;

								ii = PickAnOptionL (NcNaux1, Probs);

								if (ii >= L[mm].nc )		/* generate a new class */
								{
									UpdateEr = 1;
									Originallambda = L[mm].lambda [record];

									L[mm].lambda [record] = Auxvalues [ii-L[mm].nc];
									L[mm].lambdalist [ L[mm].lambdaclass[record] ] = L[mm].lambda [record];
									L[mm].nic [L[mm].lambdaclass[record]] ++;

								}
								else
								{
									if ( ii == L[mm].lambdaclass[record] )
									{
										L[mm].nic [L[mm].lambdaclass[record]] ++;
									}
									else
									{
										UpdateEr = 1;
										Originallambda = L[mm].lambda [record];

										L[mm].lambda [record] = L[mm].lambdalist[ii];
										Originalclass = L[mm].lambdaclass [record];
										L[mm].lambdaclass [record] = ii;

										free (L[mm].nic);
										L[mm].nic = (int*) calloc ( L[mm].nc-1, sizeof(int));
										for (jj=0; jj<N; jj++)
										{
											if (Y[L[mm].up].observations [jj] != Missingvalue)
											{
												if (L[mm].lambdaclass [jj] > Originalclass)
													L[mm].lambdaclass [jj] --;

												L[mm].nic [L[mm].lambdaclass[jj]] ++;
											}
										}
										Copydouble = (double*) malloc (sizeof(double) * L[mm].nc);
										for (cla=0; cla<L[mm].nc; cla++)
											if ( cla > Originalclass ) { Copydouble [ cla - 1] = L[mm].lambdalist [cla]; } else { Copydouble [cla] = L[mm].lambdalist [cla];}
										free (L[mm].lambdalist);
										Arraysize = sizeof (double) * ( L[mm].nc - 1 );
										L[mm].lambdalist = (double*) malloc ( Arraysize );
										memcpy ( L[mm].lambdalist, Copydouble, Arraysize);
										free (Copydouble);
										L[mm].nc --;
									}
								}
							}
							else
							{	/* in the case where record is not a singleton */

								NcNaux = L[mm].nc + Naux;
								Probs = (double*) calloc ( NcNaux, sizeof (double));

								for (cla=0; cla<Naux; cla++)
									Auxvalues [cla] = RNormal ( L[mm].meanlambda, sqrt ( L[mm].lambdasigma ) );
							
								for (trait=0; trait<K; trait++)
									Eri[trait] =  Y[trait].errors[record];

								for (cla=0; cla<L[mm].nc; cla++)
								{
									Eri [L[mm].down] += (Y [L[mm].up].stobs [record] * ( L[mm].lambda[record] - L[mm].lambdalist[cla] ));
									Probs [cla] += ( LikelihoodExp (K, iR, Eri, Products) );
									Eri [L[mm].down] = Y[ L[mm].down].errors [record];
								}
								for (cla=0; cla<Naux; cla++)
								{
									Eri [L[mm].down] += (Y [L[mm].up].stobs [record] * ( L[mm].lambda[record] - Auxvalues [cla] ));
									Probs [cla+L[mm].nc] += ( LikelihoodExp (K, iR, Eri, Products) );
									Eri [L[mm].down] = Y[ L[mm].down].errors [record];
								}

								for (cla=0; cla<L[mm].nc; cla++)
									Probs [cla] += log ( (double) L[mm].nic [cla] );

								temp = log ( L[mm].alpha/(double) Naux );
								for (cla=0; cla<Naux; cla++)
									Probs [cla+L[mm].nc] += temp;

								ii = PickAnOptionL (NcNaux, Probs);

								if (ii<L[mm].nc) 
								{
									if ( ii != L[mm].lambdaclass [record] )
									{
										UpdateEr = 1;
										Originallambda = L[mm].lambda [record];								
									}
									L[mm].lambdaclass [record] = ii;
									L[mm].lambda [record] = L[mm].lambdalist[ii];
									L[mm].nic [ii] ++;
								}
								else
								{	/* generate a new class */

									UpdateEr = 1;
									Originallambda = L[mm].lambda [record];

									L[mm].lambdaclass [record] = L[mm].nc;
									L[mm].lambda [record] = Auxvalues [ii-L[mm].nc];

									Arraysize = sizeof(double) * L[mm].nc;
									Copydouble = (double*) malloc ( Arraysize );
									memcpy ( Copydouble, L[mm].lambdalist, Arraysize);
									free (L[mm].lambdalist);
									L[mm].lambdalist = (double*) malloc ( sizeof(double) * (L[mm].nc+1) );
									memcpy ( L[mm].lambdalist, Copydouble, Arraysize);
									L[mm].lambdalist[L[mm].nc] = L[mm].lambda [record];
									free (Copydouble);

									Arraysize = sizeof (int) * L[mm].nc;
									Copyint = (int*) malloc ( Arraysize );
									memcpy ( Copyint, L[mm].nic, Arraysize );
									free (L[mm].nic);
									L[mm].nic = (int*) malloc ( sizeof(int) * (L[mm].nc+1) );
									memcpy ( L[mm].nic, Copyint, Arraysize );
									L[mm].nic [L[mm].nc] = 1;
									free (Copyint);

									L[mm].nc ++;

								}
							}	/* if else (singleton or not) */
							free (Probs);

							/* Update Er */
							if (UpdateEr)
								Y[L[mm].down].errors [record] += ( Y [ L[mm].up].stobs[record] * (Originallambda - L[mm].lambda [record]));

						}	/* missing or not */
					}		/* record */
					break;	/* L[mm].update==2 */
				}			/* switch */
			}				/* mm */
		}					/* SEM */

		/* Update Lambda */
		if (SEM)
		{
			for (mm=0; mm<M; mm++)
			{
				switch (L[mm].update) {
				case 1:
					L[mm].vlm[0] = 1.0 / (L[mm].y2[0] * iR [ L[mm].down * K + L[mm].down ] + 1.0/L[mm].lambdasigma);
					L[mm].lambdahat[0] = 1.0/L[mm].lambdasigma * L[mm].meanlambda;
					for (record=0; record<N; record++)
					{
						if (Y[L[mm].up].observations [record] != Missingvalue)
						{
							temp = 0.0;
							for (trait=0; trait<K; trait++)
							{
								if(L[mm].down == trait) 
								{
									temp += ( ( Y[trait].errors[record] + L[mm].lambda[record] * Y [L[mm].up].stobs [record] ) * iR [L[mm].down * K + trait] );
								}
								else
								{
									temp += ( Y[trait].errors[record] * iR [ L[mm].down * K + trait ] );
								}
							}
							temp *= Y [ L[mm].up].stobs[record];
							L[mm].lambdahat[0] += temp;
						}
					}
					L[mm].lambdahat[0] *= L[mm].vlm[0];
					temp = RNormal ( L[mm].lambdahat[0], sqrt(L[mm].vlm[0]) );

					/* update Er */
					for (record=0; record<N; record++)
					{
						if (Y[L[mm].up].observations [record] != Missingvalue)
						{
							Y[L[mm].down].errors[record] += ( Y [ L[mm].up].stobs[record] * (L[mm].lambda [record] - temp));
							L[mm].lambda [record] = temp;
						}
					}
					if (DoSampling)
						L[mm].sampledlambda[CumNs] = temp;
					break;

				case 2:
					L[mm].lambdahat = (double*) calloc (L[mm].nc, sizeof(double));
					L[mm].vlm = (double*) calloc (L[mm].nc, sizeof(double));
					for (record=0; record<N; record++)
					{
						if (Y[L[mm].up].observations [record] != Missingvalue)
						{
							cla = L[mm].lambdaclass[record];
							L[mm].vlm [cla] += L[mm].y2[record] * iR [ L[mm].down * K + L[mm].down ];
							temp=0.0;
							for (trait=0; trait<K; trait++)
							{
								if(L[mm].down == trait) 
								{
									temp += ( ( Y[trait].errors[record] + L[mm].lambda[record] * Y [L[mm].up].stobs [record] ) * iR [ L[mm].down * K + trait ] );
								}
								else
								{
									temp += ( Y[trait].errors[record] * iR [ L[mm].down * K + trait ] );
								}
							}
							temp *= Y [ L[mm].up].stobs [record];
							L[mm].lambdahat [cla] += temp;
						}
					}
					for (cla=0; cla<L[mm].nc; cla++)
					{
						L[mm].vlm [cla] = 1.0/(L[mm].vlm [cla] + 1.0/L[mm].lambdasigma);
						L[mm].lambdahat [cla] = L[mm].vlm [cla] * ( L[mm].lambdahat [cla] + L[mm].meanlambda/L[mm].lambdasigma );
						L[mm].lambdalist [cla] = RNormal ( L[mm].lambdahat[cla], sqrt( L[mm].vlm [cla]) );
					}

					/* update Er */
					for (record=0; record<N; record++)
					{
						if (Y[L[mm].up].observations [record] != Missingvalue)
						{
							temp = L[mm].lambdalist [L[mm].lambdaclass[record]];
							Y[L[mm].down].errors [record] += ( Y [ L[mm].up].stobs [record] * (L[mm].lambda [record] - temp));
							L[mm].lambda [record] = temp;

							if (DoSampling)
							{
								L[mm].sampledclass [CumNs*N+record] = L[mm].lambdaclass[record];
								L[mm].sampledlambda[CumNs*N+record] = L[mm].lambda[record];
							}
						}
					}
					free (L[mm].lambdahat); free(L[mm].vlm);
					break;

				case 3:
					temp = sqrt(BecomeVague*L[mm].sigmaPi);
					temp2 = sqrt(L[mm].sigmaPi);
					for(ii=0, sumPm2=0.0; ii<L[mm].Nsp; ii++)
					{
						MHprob=0.0;
						propPm = RNormal(L[mm].P[ii],SdPm);
						if(ii<UseVaguePrior)
						{
							MHprob += LogDNormal(propPm,0.0,temp);
							MHprob -= LogDNormal(L[mm].P[ii],0.0,temp);
						}
						else
						{
							if(SplinePrior==1)
							{
								MHprob += LogDNormal(propPm,L[mm].P[ii-1],temp2);
								MHprob -= LogDNormal(L[mm].P[ii],L[mm].P[ii-1],temp2);
							}
							else
							{
								temp = 2.0*L[mm].P[ii-1]-L[mm].P[ii-2];
								MHprob += LogDNormal(propPm,temp,temp2);
								MHprob -= LogDNormal(L[mm].P[ii],temp,temp2);
							}
						}

						for(record=0; record<L[mm].NhaveithSP[ii]; record++)
						{
							target = L[mm].HaveithSP[ii*N+record];
							for(jj=0; jj<K; jj++)
							{
								Eri[jj] = Y[jj].errors[target];
							}
							MHprob -= LikelihoodExp(K, iR, Eri, Products);
							prop = L[mm].lambda[target] + (propPm - L[mm].P[ii])*L[mm].B[target*L[mm].Nsp+ii]; /* new lambda value for targer */
							Eri[L[mm].down] += Y[L[mm].up].stobs[target] * (L[mm].lambda[target] - prop);
							MHprob += LikelihoodExp(K, iR, Eri, Products);
						}

						if(MHprob > log(rnd()))
						{
							/* update lambda and residual */
							for(record=0; record<L[mm].NhaveithSP[ii]; record++)
							{
								target = L[mm].HaveithSP[ii*N+record];
								prop = L[mm].lambda[target] + (propPm - L[mm].P[ii])*L[mm].B[target*L[mm].Nsp+ii];
								Y[L[mm].down].errors[target] += Y[L[mm].up].stobs[target] * (L[mm].lambda[target] - prop);
								L[mm].lambda[target] = prop;
							}
							/* update P */
							L[mm].P[ii] = propPm;
							/* count */
							L[mm].Naccepted++;
						}

						if(ii<UseVaguePrior)
						{
							sumPm2 += pow(L[mm].P[ii],2.0)/BecomeVague;
						}
						else
						{
							if(SplinePrior==1)
							{
								sumPm2 += pow(L[mm].P[ii]-L[mm].P[ii-1],2.0);
							}
							else
							{
								sumPm2 += pow(L[mm].P[ii]-2*L[mm].P[ii-1]+L[mm].P[ii-2],2.0);
							}
						}

						if(DoSampling) L[mm].sampledPi[CumNs*L[mm].Nsp+ii] = L[mm].P[ii];
					} /* ii */
					break;

				}	/* switch*/		
			}		/* mm */			
		}			/* SEM */

		/* update of sigmaPi */
		if(SEM)
			for (mm=0; mm<M; mm++)
				if(L[mm].update==3)
					L[mm].sigmaPi = Rinchi(L[mm].Nsp+L[mm].Nupi, (sumPm2+L[mm].Nupi*L[mm].S2pi)/(L[mm].Nsp+L[mm].Nupi));

		/* Sampling of iG and iR */
		if (DoSampling)
		{
			if (CovR == 0 && Polygenic == 0)
			{
				for (trait=0; trait<K; trait++ )
				{
					SamplediR [CumNs * K + trait] = iR [trait * K + trait];
				}
			}
			else
			{
				for (trait=0; trait<K; trait++ )
				{
					if (CovR == 0) SamplediR [CumNs * K + trait] = iR [trait * K + trait];
					for (ii=trait; ii<K; ii++ )
					{
						jj = CumNs * KK12 + PositionInSym (trait, ii, K);
						if (Polygenic)  SamplediG [ jj ] = iG [trait * K + ii];
						if (CovR != 0 ) SamplediR [ jj ] = iR [trait * K + ii];
					}
				}
			}
		}

		/* update errors of missign records */
		for (trait=0; trait<K; trait++)
		{
			for (record=0; record<Y[trait].nmiss; record++)
			{
				temp = 0.0;
				for (ii=0; ii<K; ii++)
					if (ii != trait)
						temp += Y[ii].errors[Y[trait].whichrec[record]] * iR[trait*K + ii];

				temp *= (-1.0/iR[trait*K + trait]);
				Y[trait].errors [ Y[trait].whichrec[record] ] = RNormal (temp, sqrt(1.0/iR[trait*K + trait]));
			}
		}

		if ( DoSampling )
		{
			for (record=0; record<N; record++)
			{
				for (ii=0; ii<K; ii++)
				{
					Eri[ii] =  Y[ii].errors[record];
					if (AddtoMean) Y[ii].meanerrors[record] += Y[ii].errors[record];
				}
				SampledLikelihood [CumNs*N+record] = LikeConstant*pow(DetR,-0.5)*exp(LikelihoodExp(K, iR, Eri, Products));
			}
		}

	}	/* mcmc repeat ends */

	putchar ('\n');
	if(CovR==1) printf ( "Acceptance rate of residual var.: %f\n", (double)NacceptedR/((double)(Totalite*K)));
	if(SEM)
		for(mm=0; mm<M; mm++)
			if(L[mm].update==3)
				printf("Acceptance rate of P: %f\n", (double)L[mm].Naccepted/((double)(Totalite*L[mm].Nsp)));

	printf ( "MCMC iterations are finished. Free memory...");

	/* Free */
	if(Polygenic)
	{
		for (additive=0; additive<PK; additive++)
			free(Z[additive].whichrec);
		free (iG); free(Saproducts); free(Sa); free(iSa);
	}

	if (F>0) 
		for (ii=0; ii<F; ii++) 
			free (X[ii].whichrec);

	if (SEM)
	{
		for (mm=0; mm<M; mm++)
		{
			free(L[mm].lambda); 
			switch (L[mm].update) {
			case 1:
				free(L[mm].lambdahat); free(L[mm].vlm);free(L[mm].y2);
				break;
			case 2:
				free(L[mm].lambdaclass); free(L[mm].lambdalist);free(L[mm].nic); free (Auxvalues);free(L[mm].y2);
				break;
			case 3:
				free(L[mm].HaveithSP); free(L[mm].NhaveithSP); free(L[mm].P); free(L[mm].knot);
				break;
			}
		}	
	}

	for (trait=0; trait<K; trait++) { free(Y[trait].errors); free(Y[trait].stobs); if(Y[trait].nmiss>0) free(Y[trait].whichrec);}

	free(Dummy); free(Products); free(Eri); 
	free(iR); free(R); free(Se); free(iSe);
	if(CovR==1){free(propR);free(propiR);}

	printf ( "done\n");
	return(0);
}
